Náš dataset obsahuje informácie o ovzduší a meracích staniciach ktoré
dané merania vykonávali. Je zložený z dvoch súborov a to: -
measurements.csv - stations.csv ## Skratky: - PM2.5 - Particulate Matter
(µg/m3)
- PM10 - Particulate Matter (µg/m3)
- NOx - Nitrogen
Oxides (µg/m3)
- NO2 - Nitrogen Dioxide (µg/m3)
- SO2 - Sulfur
Dioxide (µg/m3)
- CO - Carbon Monoxide emissions (µg/m3)
- CO2
- Carbon Dioxide (µg/m3)
- PAHs - Polycyclic Aromatic Hydrocarbons
(µg/m3)
- NH3 - Ammonia trace (µg/m3)
- Pb - Lead (µg/m3)
- TEMP - Temperature (degree Celsius)
- DEWP - Dew point
temperature (degree Celsius)
- PRES - Pressure (hPa, <100,
1050>)
- RAIN - Rain (mm)
- WSPM - Wind Speed (m/s)
-
WD - Wind Direction
- VOC - Volatile Organic Compounds
- CFCs
- Chlorofluorocarbons
- C2H3NO5 - Peroxyacetyl nitrate
- H2CO
- Plywood emit formaldehyde
- GSTM1 - Glutathione-S transferase M1
- 1-OHP - 1-hydroxypyrene
- 2-OHF - 2-hydroxyfluorene
-
2-OHNa - 2-hydroxynaphthalene
- N2 - Nitrogen
- O2 - Oxygen
- O3 - Ozone
- Ar - Argon
- Ne - Neon
- CH4 -
Methane
- He - Helium
- Kr - Krypton
- I2 - Iodine
-
H2 - Hydrogen
- Xe - Xenon
getwd()
## [1] "/Users/roman/Desktop/development/oznal_project1/part2"
#setwd('/Users/roman.osadsky/Desktop/ROMAN/FIIT/OZNAL_PROJECT/oznal_project1')
setwd('/Users/roman/Desktop/development/oznal_project1/')
#setwd(getwd())
Importovanie knižníc
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ✔ readr 2.1.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(caret)
## Loading required package: lattice
##
## Attaching package: 'caret'
##
## The following object is masked from 'package:purrr':
##
## lift
library(tidyr)
library(e1071)
library(corrplot)
## corrplot 0.92 loaded
library(MASS)
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
library(reshape2)
##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smiths
library(gridExtra)
##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
##
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(glmnet)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loaded glmnet 4.1-8
library(PRROC)
Načítanie dát
labor_measurements <- read.csv('/Users/roman/Desktop/development/oznal_project1/data/measurements.csv', sep = '\t')
labor_stations <- read.csv('/Users/roman/Desktop/development/oznal_project1/data/stations.csv', sep = '\t')
Zobrazenie prvých riadkov datasetov, ich rozmerov a typov atribútov
head(labor_measurements)
## latitude PM2.5 NOx PM10 C2H3NO5 CH4 longitude Pb
## 1 37.71715 8.47714 9.21522 9.38738 1.51791 7.84989 -122.40433 59.51096
## 2 -1.26753 8.95287 8.04229 8.78120 0.86506 7.48764 116.82887 57.03798
## 3 25.87498 11.41105 7.80172 6.67488 0.20478 8.32353 86.59611 78.98660
## 4 53.02330 5.41284 9.53800 7.01510 0.13092 6.36073 -1.48119 60.70261
## 5 -7.09810 5.59043 6.87675 6.19817 0.19487 6.55752 109.32430 79.16126
## 6 -0.35817 6.91056 6.76613 10.57920 4.65631 6.57000 42.54536 54.80587
## TEMP NH3 SO2 O3 CO PAHs H2CO PRES warning
## 1 20.05101 10.43604 5.81201 7.77502 9.69678 8.62090 47.64810 1139.127 0
## 2 20.34841 9.19522 8.62588 11.61089 4.92266 6.96953 42.66369 1144.575 1
## 3 19.57588 5.00666 6.37696 11.20430 7.61260 7.88382 48.15920 1131.932 1
## 4 20.51904 8.38101 9.96460 6.31462 6.68230 9.13018 30.33567 1176.702 0
## 5 -1.60515 8.05129 9.76426 6.23701 4.50154 8.72864 49.04137 1058.874 0
## 6 3.27832 9.64108 9.68614 11.31083 7.97488 8.54959 34.31598 1078.051 1
## CFCs
## 1 74.87342
## 2 70.13254
## 3 70.76611
## 4 68.30109
## 5 70.50372
## 6 81.84559
Počet riadkov: labor_measurements
nrow(labor_measurements)
## [1] 12118
Počet stĺpcov: labor_measurements
ncol(labor_measurements)
## [1] 18
str(labor_measurements)
## 'data.frame': 12118 obs. of 18 variables:
## $ latitude : num 37.72 -1.27 25.87 53.02 -7.1 ...
## $ PM2.5 : num 8.48 8.95 11.41 5.41 5.59 ...
## $ NOx : num 9.22 8.04 7.8 9.54 6.88 ...
## $ PM10 : num 9.39 8.78 6.67 7.02 6.2 ...
## $ C2H3NO5 : num 1.518 0.865 0.205 0.131 0.195 ...
## $ CH4 : num 7.85 7.49 8.32 6.36 6.56 ...
## $ longitude: num -122.4 116.83 86.6 -1.48 109.32 ...
## $ Pb : num 59.5 57 79 60.7 79.2 ...
## $ TEMP : num 20.05 20.35 19.58 20.52 -1.61 ...
## $ NH3 : num 10.44 9.2 5.01 8.38 8.05 ...
## $ SO2 : num 5.81 8.63 6.38 9.96 9.76 ...
## $ O3 : num 7.78 11.61 11.2 6.31 6.24 ...
## $ CO : num 9.7 4.92 7.61 6.68 4.5 ...
## $ PAHs : num 8.62 6.97 7.88 9.13 8.73 ...
## $ H2CO : num 47.6 42.7 48.2 30.3 49 ...
## $ PRES : num 1139 1145 1132 1177 1059 ...
## $ warning : num 0 1 1 0 0 1 0 0 0 1 ...
## $ CFCs : num 74.9 70.1 70.8 68.3 70.5 ...
head(labor_stations)
## latitude revision longitude code QoS location
## 1 -1.60000 2019-02-02 103.61667 ID maitennce Asia/Jakarta
## 2 29.84576 06/06/2022, 00:00:00 -90.10674 US average America/Chicago
## 3 33.03699 24 Feb 2015 -117.29198 US good America/Los_Angeles
## 4 31.84568 2017/03/10 -102.36764 US accep America/Chicago
## 5 16.06213 2015/05/29 76.05860 IN average Asia/Kolkata
## 6 19.04222 04/18/2013, 00:00:00 -98.11889 MX accep America/Mexico_City
Počet riadkov: labor_stations
nrow(labor_stations)
## [1] 1061
Počet stĺpcov: labor_stations
ncol(labor_stations)
## [1] 6
str(labor_stations)
## 'data.frame': 1061 obs. of 6 variables:
## $ latitude : num -1.6 29.8 33 31.8 16.1 ...
## $ revision : chr "2019-02-02" "06/06/2022, 00:00:00" "24 Feb 2015" "2017/03/10" ...
## $ longitude: num 103.6 -90.1 -117.3 -102.4 76.1 ...
## $ code : chr "ID" "US" "US" "US" ...
## $ QoS : chr "maitennce" "average" "good" "accep" ...
## $ location : chr "Asia/Jakarta" "America/Chicago" "America/Los_Angeles" "America/Chicago" ...
Ked sa pozrieme na labor_stations tam už to tak nieje preto si skontrolujeme formáty
# Kontrola jedinečných hodnôt a ich počtu
qos_count <- table(labor_stations$QoS)
#qos_count
# Nahradenie špecifických hodnôt
labor_stations$QoS[labor_stations$QoS == "acceptable"] <- "accep"
labor_stations$QoS[labor_stations$QoS == "maitennce"] <- "maintenance"
# Opätovná kontrola jedinečných hodnôt a ich počtov po nahradení
qos_count <- table(labor_stations$QoS)
#qos_count
helper
parse_dates <- function(date_string) {
formats <- c("%Y/%m/%d", "%d %b %Y", "%Y-%m-%d", "%m/%d/%Y, %H:%M:%S", "%d/%m/%Y, %H:%M:%S")
for (format in formats) {
parsed_date <- as.POSIXct(date_string, format = format, tz = "")
if (!is.na(parsed_date)) {
# urob datum na format "28 Aug 2016"
return(format(parsed_date, "%d %b %Y"))
}
}
return(NA)
}
# Kontrola jedinečných hodnôt a ich počtu
revision_count <- table(labor_stations$revision)
#revision_count
# Nahradenie špecifických hodnôt
labor_stations$revision <- sapply(labor_stations$revision, parse_dates)
cat('Identifikované rozdielne formáty dátumu sme zjednotili na formát dd mmm yyyy')
## Identifikované rozdielne formáty dátumu sme zjednotili na formát dd mmm yyyy
# Opätovná kontrola jedinečných hodnôt a ich počtov po nahradení
revision_count <- table(labor_stations$revision)
#revision_count
# Nahradenie prázdnych reťazcov reťazcami NA
labor_measurements[labor_measurements == ""] <- NA
labor_stations[labor_stations == ""] <- NA
Náhľad nato koľko dát chýba
missing_data_percentage <- sum(is.na(labor_measurements)) / nrow(labor_measurements) * 100
cat(sprintf("Chýbajúce dáta v labor_measurements tvoria %.5f%% dát", missing_data_percentage))
## Chýbajúce dáta v labor_measurements tvoria 6.03235% dát
cat("\n")
missing_data_percentage_stations <- sum(is.na(labor_stations)) / nrow(labor_stations) * 100
cat(sprintf("Chýbajúce dáta v labor_stations tvoria %.5f%% dát", missing_data_percentage_stations))
## Chýbajúce dáta v labor_stations tvoria 0.09425% dát
cat("\n")
Spočítanie a vypísanie riadkov s hodnotami NA pre labor_measurements
na_count_measurements <- sum(!complete.cases(labor_measurements))
cat('NaN hodnoty pre labor_measurements:', na_count_measurements, '\n')
## NaN hodnoty pre labor_measurements: 712
Spočítanie a vypísanie riadkov s hodnotami NA pre labor_stations
na_count_stations <- sum(!complete.cases(labor_stations))
cat('NaN hodnoty pre labor_stations:', na_count_stations, '\n')
## NaN hodnoty pre labor_stations: 1
Prázdne hodnoty sme sa rozhodli dropnúť
labor_measurements <- na.omit(labor_measurements)
labor_stations <- na.omit(labor_stations)
if (any(duplicated(labor_measurements))) {
duplicates_count <- sum(duplicated(labor_measurements))
cat('labor_measurements pocet duplikatov', duplicates_count, '\n')
} else {
cat('labor_measurements no duplicates\n')
}
## labor_measurements pocet duplikatov 168
if (any(duplicated(labor_stations))) {
duplicates_count_stations <- sum(duplicated(labor_stations))
cat('labor_stations pocet duplikatov', duplicates_count_stations, '\n')
} else {
cat('labor_stations no duplicates\n')
}
## labor_stations no duplicates
Dropnutie duplicitných záznamov
labor_measurements <- unique(labor_measurements)
labor_stations <- unique(labor_stations)
# Odstránenie nepotrebných stĺpcov
labor_stations <- labor_stations[, -4]
#head(labor_stations)
# Merge
df <- merge(labor_measurements, labor_stations, by = c("latitude", "longitude"), all = FALSE)
# Odstránenie stĺpcov 'latitude' a 'longitude' lebo ich už viac nepotrebujeme
df <- df[, !names(df) %in% c('latitude', 'longitude')]
# reorganizacia stlpcov
df <- df[, c('location', 'warning', 'QoS', 'revision', 'TEMP', 'PRES', 'PM2.5', 'NOx', 'PM10', 'C2H3NO5', 'CH4', 'Pb', 'NH3', 'SO2', 'O3', 'CO', 'PAHs', 'H2CO', 'CFCs')]
head(df)
## location warning QoS revision TEMP PRES PM2.5
## 1 Africa/Mogadishu 0 accep 27 Jan 2018 39.88329 1154.420 5.72081
## 2 Africa/Mogadishu 0 maintenance 01 Nov 2018 39.88329 1154.420 5.72081
## 3 Africa/Mogadishu 0 average 12 Nov 2014 39.88329 1154.420 5.72081
## 4 Africa/Mogadishu 0 excellent 26 Dec 2019 39.88329 1154.420 5.72081
## 5 Africa/Mogadishu 0 accep 27 Jan 2018 15.78476 1191.586 6.15657
## 6 Africa/Mogadishu 0 maintenance 01 Nov 2018 15.78476 1191.586 6.15657
## NOx PM10 C2H3NO5 CH4 Pb NH3 SO2 O3 CO
## 1 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212 11.00024 8.05639 7.00959
## 2 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212 11.00024 8.05639 7.00959
## 3 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212 11.00024 8.05639 7.00959
## 4 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212 11.00024 8.05639 7.00959
## 5 7.90747 7.83837 0.46016 7.94379 40.92438 4.90186 9.26682 6.62422 5.60047
## 6 7.90747 7.83837 0.46016 7.94379 40.92438 4.90186 9.26682 6.62422 5.60047
## PAHs H2CO CFCs
## 1 9.05202 34.55293 68.66001
## 2 9.05202 34.55293 68.66001
## 3 9.05202 34.55293 68.66001
## 4 9.05202 34.55293 68.66001
## 5 6.87025 54.58942 68.03564
## 6 6.87025 54.58942 68.03564
df_out <- df[, c('TEMP', 'PRES', 'PM2.5', 'NOx', 'PM10', 'C2H3NO5', 'CH4', 'Pb', 'NH3', 'SO2', 'O3', 'CO', 'PAHs', 'H2CO', 'CFCs')]
head(df_out)
## TEMP PRES PM2.5 NOx PM10 C2H3NO5 CH4 Pb NH3
## 1 39.88329 1154.420 5.72081 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212
## 2 39.88329 1154.420 5.72081 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212
## 3 39.88329 1154.420 5.72081 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212
## 4 39.88329 1154.420 5.72081 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212
## 5 15.78476 1191.586 6.15657 7.90747 7.83837 0.46016 7.94379 40.92438 4.90186
## 6 15.78476 1191.586 6.15657 7.90747 7.83837 0.46016 7.94379 40.92438 4.90186
## SO2 O3 CO PAHs H2CO CFCs
## 1 11.00024 8.05639 7.00959 9.05202 34.55293 68.66001
## 2 11.00024 8.05639 7.00959 9.05202 34.55293 68.66001
## 3 11.00024 8.05639 7.00959 9.05202 34.55293 68.66001
## 4 11.00024 8.05639 7.00959 9.05202 34.55293 68.66001
## 5 9.26682 6.62422 5.60047 6.87025 54.58942 68.03564
## 6 9.26682 6.62422 5.60047 6.87025 54.58942 68.03564
outliers_list <- lapply(df_out, function(column) {
IQR_value <- IQR(column, na.rm = TRUE)
upper_bound <- quantile(column, 0.75, na.rm = TRUE) + 1.5 * IQR_value
lower_bound <- quantile(column, 0.25, na.rm = TRUE) - 1.5 * IQR_value
lower_outliers <- which(column < lower_bound)
upper_outliers <- which(column > upper_bound)
return(list(lower_outliers = lower_outliers, upper_outliers = upper_outliers))
})
plot(df_out[[1]], main=paste("Outliers in", names(df_out)[1]), col="blue", pch=19)
if(length(outliers_list[[1]]$lower_outliers) > 0) {
points(outliers_list[[1]]$lower_outliers, df_out[outliers_list[[1]]$lower_outliers, 1], col="red", pch=19)
}
if(length(outliers_list[[1]]$upper_outliers) > 0) {
points(outliers_list[[1]]$upper_outliers, df_out[outliers_list[[1]]$upper_outliers, 1], col="red", pch=19)
}
numeric_df <- df[sapply(df, is.numeric)]
ggpairs(numeric_df,
lower = list(continuous = "smooth"),
diag = list(continuous = "barDiag"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
set.seed(123)
sample <- sample(c(TRUE, FALSE), nrow(df), replace = TRUE, prob = c(0.5, 0.3))
train_data <- df[sample,]
test_data <- df[!sample,]
# PM2.5 – the most harmful pollution
# PM10 - Particulate Matter (µg/m3)
# SO2 - Sulfur Dioxide (µg/m3)
# CO - Carbon Monoxide emissions (µg/m3)
PM25Model <- lm(PM2.5 ~ PM10 + SO2 + CO, data = train_data)
summary(PM25Model)
##
## Call:
## lm(formula = PM2.5 ~ PM10 + SO2 + CO, data = train_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.1487 -1.0612 0.0896 1.1008 6.5548
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.834566 0.074879 171.41 <2e-16 ***
## PM10 -0.495697 0.007707 -64.32 <2e-16 ***
## SO2 -0.236042 0.006927 -34.08 <2e-16 ***
## CO 0.111353 0.007168 15.53 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.588 on 15103 degrees of freedom
## Multiple R-squared: 0.3233, Adjusted R-squared: 0.3231
## F-statistic: 2405 on 3 and 15103 DF, p-value: < 2.2e-16
par(mar = c(1, 1, 1, 1))
plot(PM25Model)
Residuals vs Fitted - Červená čiara by mala horizontálna najviac ako je možná - Ak je čiara nakrivená može to znamenať nelineárny vzťah - Detekuje odlahlé body
Normal QQ - Kontrolujeme či residuals majú normálovu distribúciu - Ak sú rezíduá normálne distribuované, body na grafe by mali vytvoriť približne priamu čiaru - Ak sú body na Q-Q grafe zakrivené alebo majú iné neštandardné vzory, môže to naznačovať, že rezíduá majú neobvyklé rozdelenie - V našom prípade to vyzerá že residualy majú normálovú distribúciu
Scale-Location (Spread location) - Čiara je takmer horizontálna takže mohli by sme povedať že residuals sú rovnomerne rozložené
Residuals vs Leverage: - zobrazuje ovplivnujúce časti - niečo podobné ako outliery - ak by sme videli v rohoch čiary tak je možné že to nieje v poriadku
predicted_values <- predict(PM25Model, newdata = test_data)
# PM25 = 12.88 - 0.11 * PM10 - 0.23 * SO2 + 0.10 * CO
cat("Expected: 5.72081 , Predicted: ",predict(PM25Model,data.frame(PM10=5.85838, SO2=11.00024, CO=7.00959)))
## Expected: 5.72081 , Predicted: 8.114608
mean_absolute_error <- mean(abs(predicted_values - test_data$PM2.5))
mean_squared_error <- mean((predicted_values - test_data$PM2.5)^2)
root_mean_squared_error <- sqrt(mean_squared_error)
data.frame(Metric = c("Mean Absolute Error", "Mean Squared Error", "Root Mean Squared Error"), Value = c(mean_absolute_error,mean_squared_error,root_mean_squared_error))
## Metric Value
## 1 Mean Absolute Error 1.255622
## 2 Mean Squared Error 2.448739
## 3 Root Mean Squared Error 1.564845
MAE - meria absolutnu chybu medzi predikovanými hodnotami a skutočnímy
MSE - meria priemernu velkost chyb p redikciach modelu, a tak kvantifikuje rozdiel medzi skutočnými a predikovanými hodnotami
RMSE -> poskytuje prehlad o tom ako dobre model predikuje alebo ako blízko sú predikcie k skutočním hodnotám
plot(x=predicted_values, y=test_data$PM2.5,
xlab='Predicted Values',
ylab='Actual Values',
main='Predicted vs. Actual Values')
head(data.frame(actual=test_data$PM2.5, predicted=predicted_values))
## actual predicted
## 2 5.72081 8.114608
## 4 5.72081 8.114608
## 5 6.15657 7.385383
## 8 6.15657 7.385383
## 11 5.01126 7.391529
## 13 5.82819 7.215109
Warning je kategorický atribút, preto ho musíme pretypovať na faktor
df$warning <- as.factor(df$warning)
set.seed(123)
#cat('nrow', nrow(df))
#cat('ncol', ncol(df))
index <- sample(seq_len(nrow(df)), round(0.7 * nrow(df)))
trainData <- df[index,]
testData <- df[-index,]
# binomial distribúcii, ktorá je vhodná pre binárnu logistickú regresiu.
logisticModel <- glm(warning ~ PM2.5 + PM10 + NH3 + C2H3NO5, data = trainData, family = binomial())
summary(logisticModel)
##
## Call:
## glm(formula = warning ~ PM2.5 + PM10 + NH3 + C2H3NO5, family = binomial(),
## data = trainData)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.130687 0.180718 -28.391 <2e-16 ***
## PM2.5 0.307108 0.011165 27.505 <2e-16 ***
## PM10 0.577838 0.016466 35.093 <2e-16 ***
## NH3 -0.177567 0.009522 -18.649 <2e-16 ***
## C2H3NO5 0.168006 0.019614 8.566 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 22738 on 16884 degrees of freedom
## Residual deviance: 19630 on 16880 degrees of freedom
## AIC: 19640
##
## Number of Fisher Scoring iterations: 6
# výpočet predpovedí logistickej regresie na základe predtým natrénovaného modelu
probabilities <- predict(logisticModel, newdata = testData, type = "response")
# Na základe vrátených pravdepodobností priraďujeme triedy
# Ak je pravdepodobnosť väčšia ako 0.5, priradí sa hodnota reprezentujúca TRUE warning
# Ak je menšia alebo rovná 0.5, priradí sa FALSE warning
glmPredictions <- ifelse(probabilities > 0.5, levels(testData$warning)[2], levels(testData$warning)[1])
# konvertuje predpovedané hodnoty na faktorovú premennú s úrovňami definovanými v originálnej premennej
glmPredictions <- factor(glmPredictions, levels = levels(testData$warning))
confMatrix <- table(Predicted = glmPredictions, Actual = testData$warning)
print(confMatrix)
## Actual
## Predicted 0 1
## 0 1542 604
## 1 1391 3699
confusionMatrix je tabuľka, ktorá ukazuje počty správnych a nesprávnych predikcií modelu v porovnaní s aktuálnymi hodnotami. V tomto prípade sú predictions predpovedané hodnoty z modelu, a testData$warning sú skutočné hodnoty cieľovej premennej. Objekt confMat potom obsahuje túto metriku spolu s ďalšími štatistikami výkonu modelu.
confMat <- confusionMatrix(glmPredictions, testData$warning)
confMatMatrix <- as.matrix(confMat$table)
Melt mení confMatMatrix do dlhého formátu. Dlhý formát je užitočný pre vizualizáciu, pretože každý riadok reprezentuje jednu hodnotu s príslušnými identifikátormi. Výsledný objekt confMatMelted má tri stĺpce: - jeden pre skutočné hodnoty - jeden pre predpovedané hodnoty - a jeden pre počet pozorovaní pre každú kombináciu skutočných a predpovedaných hodnôt
confMatMelted <- melt(confMatMatrix)
# Actual - predstavuje skutočné hodnoty
# Predicted - predikované hodnoty
# Value - počet prípadov pre každú kombináciu skutočnej a predpovedanej hodnoty
colnames(confMatMelted) <- c('Actual', 'Predicted', 'Value')
# Vytvorenie ROC objektu
roc_obj <- roc(testData$warning, probabilities)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
# Výpis ROC krivky
plot(roc_obj, main="ROC Krivka")
# Hľadanie optimálneho cut-off bodu
# Zvyčajne chceme maximalizovať Youdenov index (J), ktorý je definovaný ako sensitivity + specificity - 1
coords(roc_obj, "best", best.method="youden")
## threshold specificity sensitivity
## 1 0.5518007 0.6137061 0.7859633
Vráti napríklad: - threshold (cut-off bod) - je hodnota, ktorá slúži
ako hranica na rozhodovanie medzi dvoma možnosťami
- specificity -
pravdepodobnosť, že model predpovedá pozitívny výsledok pozorovania, keď
je výsledok skutočne pozitívny
- sensitivity - pravdepodobnosť, že
model predpovedá negatívny výsledok pozorovania, keď je výsledok
skutočne negatívny
# Vypočítajme krivku ROC a nájdime optimálny point
roc_obj <- roc(testData$warning, probabilities)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
optimal <- coords(roc_obj, "best", best.method="youden")
optimal_threshold <- optimal$threshold
# Kontrola platnosti optimálnych súradníc
if (!is.null(optimal$x) && !is.null(optimal$y)) {
# Plot ROC curve
plot(roc_obj, main="ROC Curve with Optimal Threshold")
# bod pre optimal threshold
points(optimal$x, optimal$y, pch=19, col="red")
# pridaj text label pre optimal threshold
text(optimal$x, optimal$y, labels=paste("Threshold:", round(optimal_threshold, 3)), pos=4)
segments(optimal$x, optimal$y, optimal$x, 0, col="red", lty=2)
} else {
cat("Optimal coordinates are not available.\n")
}
## Optimal coordinates are not available.
p <- ggplot(confMatMelted, aes(x = Actual, y = Predicted, fill = Value)) +
geom_tile(color = "white") + # Tiles for each combination of actual and predicted
geom_text(aes(label = Value), color = "black", size = 5) + # prida text pre kazdu dlaždicu
scale_fill_gradient(low = "white", high = "steelblue") + # farba gradientu
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1)) + # prida x-axis text
xlab('Predicted') + # X-axis label
ylab('Actual') + # Y-axis label
ggtitle('Confusion Matrix - Logistic Regression') # title
print(p)
## Lasso regresion
head(df)
## location warning QoS revision TEMP PRES PM2.5
## 1 Africa/Mogadishu 0 accep 27 Jan 2018 39.88329 1154.420 5.72081
## 2 Africa/Mogadishu 0 maintenance 01 Nov 2018 39.88329 1154.420 5.72081
## 3 Africa/Mogadishu 0 average 12 Nov 2014 39.88329 1154.420 5.72081
## 4 Africa/Mogadishu 0 excellent 26 Dec 2019 39.88329 1154.420 5.72081
## 5 Africa/Mogadishu 0 accep 27 Jan 2018 15.78476 1191.586 6.15657
## 6 Africa/Mogadishu 0 maintenance 01 Nov 2018 15.78476 1191.586 6.15657
## NOx PM10 C2H3NO5 CH4 Pb NH3 SO2 O3 CO
## 1 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212 11.00024 8.05639 7.00959
## 2 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212 11.00024 8.05639 7.00959
## 3 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212 11.00024 8.05639 7.00959
## 4 7.95611 5.85838 0.15631 9.75768 40.34402 8.07212 11.00024 8.05639 7.00959
## 5 7.90747 7.83837 0.46016 7.94379 40.92438 4.90186 9.26682 6.62422 5.60047
## 6 7.90747 7.83837 0.46016 7.94379 40.92438 4.90186 9.26682 6.62422 5.60047
## PAHs H2CO CFCs
## 1 9.05202 34.55293 68.66001
## 2 9.05202 34.55293 68.66001
## 3 9.05202 34.55293 68.66001
## 4 9.05202 34.55293 68.66001
## 5 6.87025 54.58942 68.03564
## 6 6.87025 54.58942 68.03564
library(glmnet)
sample2 <- sample(c(TRUE, FALSE), nrow(df), replace = TRUE, prob = c(0.5, 0.3))
train_data <- df[sample2,]
test_data <- df[!sample2,]
# Nastavenie cieľovej hodnoty na PM2.5
y_train_data <- train_data$PM2.5
y_test_data <- test_data$PM2.5
# Nastavenie nezávislých premenných na PM10, SO2 a CO a vytvorenie matice dát
x_train_data <- data.matrix(train_data[, c('PM10', 'SO2', 'CO')])
x_test_data <- data.matrix(test_data[, c('PM10', 'SO2', 'CO')])
x_train_data_SO2_CO <- data.matrix(train_data[, c('SO2', 'CO')])
x_test_data_SO2_CO <- data.matrix(test_data[, c('SO2', 'CO')])
x_train_data_C2H3NO5_NH3_PM10 <- data.matrix(train_data[, c('C2H3NO5', 'NH3', 'PM10')])
x_test_data_C2H3NO5_NH3_PM10<- data.matrix(test_data[, c('C2H3NO5', 'NH3', 'PM10')])
# Vytvorenie krížovej validácie
# cv = Cros validation pre nájdenie optimálnej hodnoty pre Lambdu
# gausiann -> dávame vedieť že robíme lineárnu regresiu
# alpha = 1 -> pre lasso regressiu
cv_model <- cv.glmnet(x_train_data, y_train_data, alpha = 1, family = "gaussian")
cv_model_SO2_CO <- cv.glmnet(x_train_data_SO2_CO, y_train_data, alpha = 1, family = "gaussian")
cv_model_C2H3NO5_NH3_PM10 <- cv.glmnet(x_train_data_C2H3NO5_NH3_PM10, y_train_data, alpha = 1, family = "gaussian")
# Určenie optimálnej hodnoty lambda z krížovej validácie
# hyperparameter lambda (λ) - vyvažuje kompromis medzi odchýlkou a rozptylom vo výsledných koeficientoch.
# labda môže byť od 0 do positive nekonečna
best_lambda <- cv_model$lambda.min
best_lambda_SO2_CO <- cv_model_SO2_CO$lambda.min
best_lambda_C2H3NO5_NH3_PM10 <- cv_model_C2H3NO5_NH3_PM10$lambda.min
best_lambda
## [1] 0.003832762
best_lambda_SO2_CO
## [1] 0.007008515
best_lambda_C2H3NO5_NH3_PM10
## [1] 0.005560681
lambda_values <- data.frame(
Model = c("Model 1", "Model SO2_CO", "Model C2H3NO5_NH3_PM10"),
Labdas = c(best_lambda, best_lambda_SO2_CO, best_lambda_C2H3NO5_NH3_PM10)
)
lambda_values
## Model Labdas
## 1 Model 1 0.003832762
## 2 Model SO2_CO 0.007008515
## 3 Model C2H3NO5_NH3_PM10 0.005560681
# Vykreslenie grafu krížovej validácie
# Na tom plote môžeme najsť najnižšiu labdu, s najnizším cross-valudation errorom
plot(cv_model)
plot(cv_model_SO2_CO)
plot(cv_model_C2H3NO5_NH3_PM10)
# Získanie koeficientov najlepšieho Lasso regresného modelu
best_model <- glmnet(x_train_data, y_train_data, alpha = 1, lambda = best_lambda)
best_model_SO2_CO <- glmnet(x_train_data_SO2_CO, y_train_data, alpha = 1, lambda = best_lambda_SO2_CO)
best_model_C2H3NO5_NH3_PM10 <- glmnet(x_train_data_C2H3NO5_NH3_PM10, y_train_data, alpha = 1, lambda = best_lambda_C2H3NO5_NH3_PM10)
# Koeficienty linarnej regressie
# PM25 = 12.88 - 0.11 * PM10 - 0.23 * SO2 + 0.10 * CO
coef(best_model)
## 4 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 12.8818570
## PM10 -0.4989048
## SO2 -0.2330829
## CO 0.1056300
coef(best_model_SO2_CO)
## 3 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 11.049899713
## SO2 -0.363095121
## CO -0.005586441
coef(best_model_C2H3NO5_NH3_PM10)
## 4 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 12.15689355
## C2H3NO5 -0.02254332
## NH3 -0.01331249
## PM10 -0.52098774
Môžeme vidieť že pri modely s prediktormi SO2, CO má C0 bola eliminovaná z modelu čé naznačuje, že CO nemá štatisticky významný vzťah.
# Predikcia hodnôt pomocou najlepšieho modelu
# do modelu dávame lambdu ako penalizačný parameter pre minimalizáciu chyby predikcie ako prevenciu nadmerného prispôsobenia sa dátam.
y_predicted <- predict(best_model, s = best_lambda, newx = x_test_data)
y_predicted_S02_CO <- predict(best_model_SO2_CO, s = best_lambda_SO2_CO, newx = x_test_data_SO2_CO)
y_predicted_C2H3NO5_NH3_PM10 <- predict(best_model_C2H3NO5_NH3_PM10, s = best_lambda_C2H3NO5_NH3_PM10, newx = x_test_data_C2H3NO5_NH3_PM10)
sst <- sum((y_test_data - mean(y_test_data))^2)
sse <- sum((y_predicted - y_test_data)^2)
rsq <- 1 - sse/sst
sst <- sum((y_test_data - mean(y_test_data))^2)
sse <- sum((y_predicted_S02_CO - y_test_data)^2)
rsq_SO2_CO <- 1 - sse/sst
sst <- sum((y_test_data - mean(y_test_data))^2)
sse <- sum((y_predicted_C2H3NO5_NH3_PM10 - y_test_data)^2)
rsq_C2H3NO5_NH3_PM10 <- 1 - sse/sst
mse <- (sum(y_predicted - y_test_data)^2/length(y_predicted))
mse_SE2_CO <- (sum(y_predicted_S02_CO - y_test_data)^2/length(y_predicted_S02_CO))
mse_C2H3NO5_NH3_PM10 <- (sum(y_predicted_C2H3NO5_NH3_PM10 - y_test_data)^2/length(y_predicted_C2H3NO5_NH3_PM10))
rsq_values <- data.frame(
Model = c("Model 1", "Model SO2_CO", "Model C2H3NO5_NH3_PM10"),
RSquared = c(rsq, rsq_SO2_CO, rsq_C2H3NO5_NH3_PM10),
MSEError = c(mse, mse_SE2_CO, mse_C2H3NO5_NH3_PM10)
)
print(rsq_values)
## Model RSquared MSEError
## 1 Model 1 0.3278211 0.10669062
## 2 Model SO2_CO 0.1331646 0.02409379
## 3 Model C2H3NO5_NH3_PM10 0.2750108 0.01198271
Ked sa pozrieme na MSE hodnotu, model C2H3NO5_NH3_PM10 funguje najlepšie z hľadiska presnosti predpovede, a to aj napriek nižšej RSQ hodnote.
Model 1 nemá najnižšiu MSE, a najvyčší rsquared môže byť vhodný model na dalšiu prácu.
# rozdelenie data setu na test a train
# nieje potrebne lebo je urobene hore
set.seed(123)
index <- sample(seq_len(nrow(df)), round(0.7 * nrow(df)))
trainData <- df[index,]
testData <- df[-index,]
# Trénovanie modelu Linear Discriminant Analysis (LDA)
ldaModel <- lda(warning ~ PM2.5 + PM10 + NH3 + C2H3NO5, data = trainData)
# Predpovedanie na testovacom súbore údajov
ldaPredictions <- predict(ldaModel, newdata = testData)$class
# evaluacia modelu
table(Predicted = ldaPredictions, Actual = testData$warning)
## Actual
## Predicted 0 1
## 0 1476 579
## 1 1457 3724
# Trénovanie modelu Linear Discriminant Analysis (LDA)
ldaModelTweak1 <- lda(warning ~ NH3 + C2H3NO5, data = trainData)
# Predpovedanie na testovacom súbore údajov
ldaPredictionsTweak1 <- predict(ldaModelTweak1, newdata = testData)$class
# evaluacia modelu
table(Predicted = ldaPredictionsTweak1, Actual = testData$warning)
## Actual
## Predicted 0 1
## 0 1 0
## 1 2932 4303
SVM sú typom algoritmu učenia pod dohľadom, ktorý sa používa na klasifikačné aj regresné úlohy. Sú obzvlášť výkonné pre vysokorozmerné údaje. Cieľom SVM je nájsť najlepšiu oddeľujúcu hyperplochu, ktorá rozdeľuje rôzne triedy v priestore príznakov. Vektory, ktoré definujú hyperplochu, sa nazývajú podporné vektory.
Kernel je metóda používaná v SVM, ktorá umožňuje učenie vo vysokodimenzionálnych priestoroch príznakov bez explicitnej transformácie údajov do týchto dimenzií. Je to užitočné najmä vtedy, keď hranica medzi triedami nie je lineárna. Použitím funkcie jadra môže SVM vykonávať zložité klasifikácie pomocou lineárneho klasifikátora implicitným mapovaním vstupných údajov do vysokodimenzionálnych priestorov príznakov.
Kernel RBF je jednou z najbežnejších funkcií používaných v SVM. Označuje sa aj ako Gaussov kernel. Výhody použitia jadra RBF v SVM - RBF jadro je veľmi flexibilné a dokáže riešiť prípady, keď je vzťah medzi značkami tried a atribútmi nelineárny. - Funguje dobre vo vysokodimenzionálnych priestoroch, aj keď je počet dimenzií väčší ako počet vzoriek. - Je účinný v scenároch, kde dátové body nie sú lineárne oddeliteľné.
kernel možnosti v e1071:
set.seed(123)
index <- sample(seq_len(nrow(df)), round(0.7 * nrow(df)))
trainData <- df[index,]
testData <- df[-index,]
# trenovanie SVM modelu
svmModel <- svm(warning ~ PM2.5 + PM10 + NH3 + C2H3NO5, data = trainData, kernel = "radial", cost = 1, scale = TRUE)
summary(svmModel)
##
## Call:
## svm(formula = warning ~ PM2.5 + PM10 + NH3 + C2H3NO5, data = trainData,
## kernel = "radial", cost = 1, scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
##
## Number of Support Vectors: 6669
##
## ( 3346 3323 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
# generovanie predpovedí
svmPredictions <- predict(svmModel, newdata = testData)
# prevod predpovedí na faktor s úrovňami ako v testovacích údajoch
svmPredictions <- factor(svmPredictions, levels = levels(testData$warning))
# Confusion Matrix pre SVM model
svmConfMatrix <- table(Predicted = svmPredictions, Actual = testData$warning)
print(svmConfMatrix)
## Actual
## Predicted 0 1
## 0 2130 400
## 1 803 3903
# pokročilé vyhodnotie
confusionMatrix(svmPredictions, testData$warning)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 2130 400
## 1 803 3903
##
## Accuracy : 0.8337
## 95% CI : (0.825, 0.8423)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6474
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.7262
## Specificity : 0.9070
## Pos Pred Value : 0.8419
## Neg Pred Value : 0.8294
## Prevalence : 0.4053
## Detection Rate : 0.2944
## Detection Prevalence : 0.3496
## Balanced Accuracy : 0.8166
##
## 'Positive' Class : 0
##
# trenovanie SVM modelu
svmModelTweak1 <- svm(warning ~ NH3 + C2H3NO5, data = trainData, kernel = "radial", cost = 1, scale = TRUE)
summary(svmModelTweak1)
##
## Call:
## svm(formula = warning ~ NH3 + C2H3NO5, data = trainData, kernel = "radial",
## cost = 1, scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
##
## Number of Support Vectors: 10379
##
## ( 5208 5171 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
# generovanie predpovedí
svmPredictionsTweak1 <- predict(svmModelTweak1, newdata = testData)
# prevod predpovedí na faktor s úrovňami ako v testovacích údajoch
svmPredictionsTweak1 <- factor(svmPredictionsTweak1, levels = levels(testData$warning))
# Confusion Matrix pre SVM model
svmConfMatrixTweak1 <- table(Predicted = svmPredictionsTweak1, Actual = testData$warning)
print(svmConfMatrixTweak1)
## Actual
## Predicted 0 1
## 0 1498 493
## 1 1435 3810
# pokročilé vyhodnotie
confusionMatrix(svmPredictionsTweak1, testData$warning)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1498 493
## 1 1435 3810
##
## Accuracy : 0.7336
## 95% CI : (0.7232, 0.7437)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4175
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.5107
## Specificity : 0.8854
## Pos Pred Value : 0.7524
## Neg Pred Value : 0.7264
## Prevalence : 0.4053
## Detection Rate : 0.2070
## Detection Prevalence : 0.2752
## Balanced Accuracy : 0.6981
##
## 'Positive' Class : 0
##
# trenovanie SVM modelu
svmModelTweak2 <- svm(warning ~ NOx + C2H3NO5 + CH4 + Pb + SO2 + O3 + CO + PAHs + H2CO + CFCs, data = trainData, kernel = "radial", cost = 1, scale = TRUE)
summary(svmModelTweak2)
##
## Call:
## svm(formula = warning ~ NOx + C2H3NO5 + CH4 + Pb + SO2 + O3 + CO +
## PAHs + H2CO + CFCs, data = trainData, kernel = "radial", cost = 1,
## scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
##
## Number of Support Vectors: 7999
##
## ( 4098 3901 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
# generovanie predpovedí
svmPredictionsTweak2 <- predict(svmModelTweak2, newdata = testData)
# prevod predpovedí na faktor s úrovňami ako v testovacích údajoch
svmPredictionsTweak2 <- factor(svmPredictionsTweak2, levels = levels(testData$warning))
# Confusion Matrix pre SVM model
svmConfMatrixTweak2 <- table(Predicted = svmPredictionsTweak2, Actual = testData$warning)
print(svmConfMatrixTweak2)
## Actual
## Predicted 0 1
## 0 2055 308
## 1 878 3995
# pokročilé vyhodnotie
confusionMatrix(svmPredictionsTweak2, testData$warning)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 2055 308
## 1 878 3995
##
## Accuracy : 0.8361
## 95% CI : (0.8274, 0.8446)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6492
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.7006
## Specificity : 0.9284
## Pos Pred Value : 0.8697
## Neg Pred Value : 0.8198
## Prevalence : 0.4053
## Detection Rate : 0.2840
## Detection Prevalence : 0.3266
## Balanced Accuracy : 0.8145
##
## 'Positive' Class : 0
##
Stupeň polynómu je kľúčovým parametrom, pretože určuje flexibilitu rozhodovacej hranice. Vyšší stupeň umožňuje, aby sa hranica viac zakrivovala, čím potenciálne vyhovuje zložitejším vzorom v údajoch. Tým sa však zvyšuje aj riziko nadmerného prispôsobenia, najmä ak je stupeň príliš vysoký vzhľadom na zložitosť alebo množstvo údajov.
# trenovanie SVM modelu
svmModelPolynomial <- svm(warning ~ PM2.5 + PM10 + NH3 + C2H3NO5, data = trainData, kernel = "polynomial", cost = 1, scale = TRUE, degree = 4)
summary(svmModelPolynomial)
##
## Call:
## svm(formula = warning ~ PM2.5 + PM10 + NH3 + C2H3NO5, data = trainData,
## kernel = "polynomial", cost = 1, degree = 4, scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 1
## degree: 4
## coef.0: 0
##
## Number of Support Vectors: 9377
##
## ( 4691 4686 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
# generovanie predpovedí
svmPredictionsPolynomial <- predict(svmModelPolynomial, newdata = testData)
# prevod predpovedí na faktor s úrovňami ako v testovacích údajoch
svmPredictionsPolynomial <- factor(svmPredictionsPolynomial, levels = levels(testData$warning))
# Confusion Matrix pre SVM model
svmConfMatrixPolynomial <- table(Predicted = svmPredictionsPolynomial, Actual = testData$warning)
print(svmConfMatrixPolynomial)
## Actual
## Predicted 0 1
## 0 1333 231
## 1 1600 4072
# pokročilé vyhodnotie
confusionMatrix(svmPredictionsPolynomial, testData$warning)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1333 231
## 1 1600 4072
##
## Accuracy : 0.747
## 95% CI : (0.7368, 0.7569)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.433
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.4545
## Specificity : 0.9463
## Pos Pred Value : 0.8523
## Neg Pred Value : 0.7179
## Prevalence : 0.4053
## Detection Rate : 0.1842
## Detection Prevalence : 0.2161
## Balanced Accuracy : 0.7004
##
## 'Positive' Class : 0
##
# trenovanie SVM modelu
svmModelLinear <- svm(warning ~ PM2.5 + PM10 + NH3 + C2H3NO5, data = trainData, kernel = "polynomial", cost = 1, scale = TRUE)
summary(svmModelLinear)
##
## Call:
## svm(formula = warning ~ PM2.5 + PM10 + NH3 + C2H3NO5, data = trainData,
## kernel = "polynomial", cost = 1, scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 1
## degree: 3
## coef.0: 0
##
## Number of Support Vectors: 9904
##
## ( 4952 4952 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
# generovanie predpovedí
svmPredictionsLinear <- predict(svmModelLinear, newdata = testData)
# prevod predpovedí na faktor s úrovňami ako v testovacích údajoch
svmPredictionsLinear <- factor(svmPredictionsLinear, levels = levels(testData$warning))
# Confusion Matrix pre SVM model
svmConfMatrixLinear <- table(Predicted = svmPredictionsLinear, Actual = testData$warning)
print(svmConfMatrixLinear)
## Actual
## Predicted 0 1
## 0 1456 190
## 1 1477 4113
# pokročilé vyhodnotie
confusionMatrix(svmPredictionsLinear, testData$warning)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1456 190
## 1 1477 4113
##
## Accuracy : 0.7696
## 95% CI : (0.7597, 0.7793)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4862
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.4964
## Specificity : 0.9558
## Pos Pred Value : 0.8846
## Neg Pred Value : 0.7358
## Prevalence : 0.4053
## Detection Rate : 0.2012
## Detection Prevalence : 0.2275
## Balanced Accuracy : 0.7261
##
## 'Positive' Class : 0
##
# Vytvorenie confusionMatrix matíc pre každý model
glmConfMatrix <- confusionMatrix(as.factor(glmPredictions), testData$warning)
ldaConfMatrix <- confusionMatrix(as.factor(ldaPredictions), testData$warning)
svmConfMatrix <- confusionMatrix(as.factor(svmPredictions), testData$warning)
ldaConfMatrixTweak1 <- confusionMatrix(as.factor(ldaPredictionsTweak1), testData$warning)
svmConfMatrixTweak1 <- confusionMatrix(as.factor(svmPredictionsTweak1), testData$warning)
svmConfMatrixTweak2 <- confusionMatrix(as.factor(svmPredictionsTweak2), testData$warning)
# Vypísanie matíc a súhrnných štatistík
print("GLM Confusion Matrix and Summary:")
## [1] "GLM Confusion Matrix and Summary:"
print(glmConfMatrix)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1542 604
## 1 1391 3699
##
## Accuracy : 0.7243
## 95% CI : (0.7138, 0.7346)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4026
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.5257
## Specificity : 0.8596
## Pos Pred Value : 0.7185
## Neg Pred Value : 0.7267
## Prevalence : 0.4053
## Detection Rate : 0.2131
## Detection Prevalence : 0.2966
## Balanced Accuracy : 0.6927
##
## 'Positive' Class : 0
##
print("LDA Confusion Matrix and Summary:")
## [1] "LDA Confusion Matrix and Summary:"
print(ldaConfMatrix)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1476 579
## 1 1457 3724
##
## Accuracy : 0.7186
## 95% CI : (0.7081, 0.729)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.3871
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.5032
## Specificity : 0.8654
## Pos Pred Value : 0.7182
## Neg Pred Value : 0.7188
## Prevalence : 0.4053
## Detection Rate : 0.2040
## Detection Prevalence : 0.2840
## Balanced Accuracy : 0.6843
##
## 'Positive' Class : 0
##
print("SVM Confusion Matrix and Summary:")
## [1] "SVM Confusion Matrix and Summary:"
print(svmConfMatrix)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 2130 400
## 1 803 3903
##
## Accuracy : 0.8337
## 95% CI : (0.825, 0.8423)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6474
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.7262
## Specificity : 0.9070
## Pos Pred Value : 0.8419
## Neg Pred Value : 0.8294
## Prevalence : 0.4053
## Detection Rate : 0.2944
## Detection Prevalence : 0.3496
## Balanced Accuracy : 0.8166
##
## 'Positive' Class : 0
##
print("LDA Tweak1 Confusion Matrix and Summary:")
## [1] "LDA Tweak1 Confusion Matrix and Summary:"
print(ldaConfMatrixTweak1)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1 0
## 1 2932 4303
##
## Accuracy : 0.5948
## 95% CI : (0.5834, 0.6061)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : 0.4955
##
## Kappa : 4e-04
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.0003409
## Specificity : 1.0000000
## Pos Pred Value : 1.0000000
## Neg Pred Value : 0.5947478
## Prevalence : 0.4053344
## Detection Rate : 0.0001382
## Detection Prevalence : 0.0001382
## Balanced Accuracy : 0.5001705
##
## 'Positive' Class : 0
##
print("SVM Tweak1 Confusion Matrix and Summary:")
## [1] "SVM Tweak1 Confusion Matrix and Summary:"
print(svmConfMatrixTweak1)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1498 493
## 1 1435 3810
##
## Accuracy : 0.7336
## 95% CI : (0.7232, 0.7437)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4175
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.5107
## Specificity : 0.8854
## Pos Pred Value : 0.7524
## Neg Pred Value : 0.7264
## Prevalence : 0.4053
## Detection Rate : 0.2070
## Detection Prevalence : 0.2752
## Balanced Accuracy : 0.6981
##
## 'Positive' Class : 0
##
print("SVM Tweak2 Confusion Matrix and Summary:")
## [1] "SVM Tweak2 Confusion Matrix and Summary:"
print(svmConfMatrixTweak2)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 2055 308
## 1 878 3995
##
## Accuracy : 0.8361
## 95% CI : (0.8274, 0.8446)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6492
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.7006
## Specificity : 0.9284
## Pos Pred Value : 0.8697
## Neg Pred Value : 0.8198
## Prevalence : 0.4053
## Detection Rate : 0.2840
## Detection Prevalence : 0.3266
## Balanced Accuracy : 0.8145
##
## 'Positive' Class : 0
##
# Výpočet ROC kriviek pre každý model
glm_roc <- roc(response = testData$warning, predictor = as.numeric(glmPredictions))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
lda_roc <- roc(response = testData$warning, predictor = as.numeric(ldaPredictions))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
svm_roc <- roc(response = testData$warning, predictor = as.numeric(svmPredictions))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
lda_rocTweak1 <- roc(response = testData$warning, predictor = as.numeric(ldaPredictionsTweak1))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
svm_rocTweak1 <- roc(response = testData$warning, predictor = as.numeric(svmPredictionsTweak1))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
svm_rocTweak2 <- roc(response = testData$warning, predictor = as.numeric(svmPredictionsTweak2))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
# Vykreslenie ROC kriviek
plot(glm_roc, main="ROC Curves", col="red")
plot(lda_roc, add=TRUE, col="blue")
plot(svm_roc, add=TRUE, col="green")
plot(lda_rocTweak1, add=TRUE, col="black")
plot(svm_rocTweak1, add=TRUE, col="pink")
plot(svm_rocTweak2, add=TRUE, col="orange")
legend("bottomright", legend=c("GLM", "LDA", "SVM", "LDA_TWEAK1", "SVM_TWEAK1", "SVM_TWEAK2"), col=c("red", "blue", "green", "black", "pink", "orange"), lwd=2)
# Pomocná funkcia na premenu zmätkovej matice na data frame vhodný pre ggplot
conf_matrix_to_df <- function(conf_matrix) {
df <- as.data.frame(as.table(conf_matrix))
colnames(df) <- c("Prediction", "Reference", "Frequency") # Meníme názvy stĺpcov pre zrozumiteľnosť a správne mapovanie
return(df)
}
# Prevod zmätkových matíc na data frame
glm_df <- conf_matrix_to_df(glmConfMatrix$table)
lda_df <- conf_matrix_to_df(ldaConfMatrix$table)
svm_df <- conf_matrix_to_df(svmConfMatrix$table)
lda_df_Tweak1 <- conf_matrix_to_df(ldaConfMatrixTweak1$table)
svm_df_Tweak1 <- conf_matrix_to_df(svmConfMatrixTweak1$table)
svm_df_Tweak2 <- conf_matrix_to_df(svmConfMatrixTweak2$table)
# Funkcia na vykreslenie heatmapy
plot_conf_matrix <- function(data, title) {
ggplot(data, aes(x = Prediction, y = Reference, fill = Frequency)) +
geom_tile(color = "white") +
scale_fill_gradient(low = "blue", high = "red") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(title = title, x = 'Predicted Class', y = 'Actual Class') +
geom_text(aes(label = Frequency), vjust = 1.5, color = "white")
}
# Vykreslenie heatmap pre každú zmätkovú maticu
plot_conf_matrix(glm_df, "GLM Confusion Matrix")
plot_conf_matrix(lda_df, "LDA Confusion Matrix")
plot_conf_matrix(svm_df, "SVM Confusion Matrix")
plot_conf_matrix(lda_df_Tweak1, "LDA Tweak1 Confusion Matrix")
plot_conf_matrix(svm_df_Tweak1, "SVM Tweak1 Confusion Matrix")
plot_conf_matrix(svm_df_Tweak2, "SVM Tweak2 Confusion Matrix")
Plocha pod krivkou PR (AUPRC) kvantifikuje celkový výkon klasifikačného modelu. Vypočítava plochu pod krivkou PR, ktorá predstavuje schopnosť modelu dosiahnuť vysokú presnosť a vysokú odvolávku súčasne pri všetkých možných prahových hodnotách. vyššia hodnota AUPRC naznačuje lepšiu výkonnosť modelu, pretože znamená väčšiu plochu pod krivkou PR, čo naznačuje, že model môže dosiahnuť vysokú presnosť a vysokú odvolávku pri rôznych prahových hodnotách. AUPRC je obzvlášť užitočná pri práci s nevyváženými súbormi údajov, kde jedna trieda prevláda oveľa viac ako druhá, pretože sa zameriava na výkonnosť pozitívnej triedy, ktorá môže byť zaujímavejšia. v súhrne krivka Precision-Recall poskytuje prehľad o kompromise medzi presnosťou a odvolávkou pri rôznych klasifikačných prahoch, zatiaľ čo AUPRC sumarizuje celkovú výkonnosť modelu pri všetkých možných prahoch. Oba sú cennými nástrojmi na hodnotenie účinnosti klasifikačných modelov, najmä v scenároch, kde je prítomná nerovnováha tried.
# Generovanie kriviek PR s parametrom krivky nastaveným na TRUE
# curve parameter musi byt true ak chem plotovat krivku
pr_glm <- pr.curve(scores.class0 = glmPredictions, weights.class0 = testData$warning == 1, curve = TRUE)
pr_lda <- pr.curve(scores.class0 = ldaPredictions, weights.class0 = testData$warning == 1, curve = TRUE)
pr_svm <- pr.curve(scores.class0 = svmPredictions, weights.class0 = testData$warning == 1, curve = TRUE)
pr_lda_Tweak1 <- pr.curve(scores.class0 = ldaPredictionsTweak1, weights.class0 = testData$warning == 1, curve = TRUE)
pr_svm_Tweak1 <- pr.curve(scores.class0 = svmPredictionsTweak1, weights.class0 = testData$warning == 1, curve = TRUE)
pr_svm_Tweak2 <- pr.curve(scores.class0 = svmPredictionsTweak2, weights.class0 = testData$warning == 1, curve = TRUE)
# kontrola štruktúri PR curve objects
str(pr_glm)
## List of 4
## $ type : chr "PR"
## $ auc.integral : num 0.716
## $ auc.davis.goadrich: num 0.716
## $ curve : num [1:609, 1:3] 1 1 1 1 1 ...
## - attr(*, "class")= chr "PRROC"
str(pr_lda)
## List of 4
## $ type : chr "PR"
## $ auc.integral : num 0.71
## $ auc.davis.goadrich: num 0.71
## $ curve : num [1:584, 1:3] 1 1 1 1 1 ...
## - attr(*, "class")= chr "PRROC"
str(pr_svm)
## List of 4
## $ type : chr "PR"
## $ auc.integral : num 0.817
## $ auc.davis.goadrich: num 0.817
## $ curve : num [1:405, 1:3] 1 1 1 1 1 ...
## - attr(*, "class")= chr "PRROC"
str(pr_lda_Tweak1)
## List of 4
## $ type : chr "PR"
## $ auc.integral : num 0.595
## $ auc.davis.goadrich: num 0.595
## $ curve : num [1:4, 1:3] 1 1 1 0 0.595 ...
## - attr(*, "class")= chr "PRROC"
str(pr_svm_Tweak1)
## List of 4
## $ type : chr "PR"
## $ auc.integral : num 0.718
## $ auc.davis.goadrich: num 0.718
## $ curve : num [1:498, 1:3] 1 1 1 1 1 ...
## - attr(*, "class")= chr "PRROC"
str(pr_svm_Tweak2)
## List of 4
## $ type : chr "PR"
## $ auc.integral : num 0.811
## $ auc.davis.goadrich: num 0.811
## $ curve : num [1:313, 1:3] 1 1 1 1 1 ...
## - attr(*, "class")= chr "PRROC"
plot(pr_glm, col = "red")
lines(pr_lda$recalls, pr_lda$precisions, col = "blue")
lines(pr_svm$recalls, pr_svm$precisions, col = "green")
legend("bottomright", legend = c("GLM", "LDA", "SVM"), col = c("red", "blue", "green"), lty = 1)
meria presnosť pozitívnych predpovedí modelu. Odpovedá na otázku: “Koľko zo všetkých prípadov, ktoré model predpovedal ako pozitívne, bolo skutočne pozitívnych? Matematicky sa presnosť vypočíta ako: Vysoká presnosť znamená, že model má tendenciu robiť menej falošne pozitívnych predpovedí.
meria schopnosť modelu správne identifikovať pozitívne prípady z celkového počtu skutočných pozitívnych prípadov. Odpovedá na otázku: “Koľko zo všetkých skutočných pozitívnych prípadov model správne predpovedal ako pozitívne?” Vysoká spätná väzba znamená, že model zachytáva väčšinu skutočných pozitívnych prípadov.
je harmonický priemer Precision a Recall, ktorý poskytuje jedinú metriku na vyváženie Precision a Recall. Je užitočné najmä vtedy, keď potrebujete rovnako zohľadniť falošne pozitívne aj falošne negatívne výsledky. SkóreF1 dosahuje najlepšiu hodnotu pri 1 (dokonalá presnosť a odvolanie) a najhoršiu pri 0.
# extrahujem precision a recall z confusion matrices
precision_glm <- glmConfMatrix$byClass['Pos Pred Value']
recall_glm <- glmConfMatrix$byClass['Sensitivity']
f1_glm <- 2 * (precision_glm * recall_glm) / (precision_glm + recall_glm)
precision_lda <- ldaConfMatrix$byClass['Pos Pred Value']
recall_lda <- ldaConfMatrix$byClass['Sensitivity']
f1_lda <- 2 * (precision_lda * recall_lda) / (precision_lda + recall_lda)
precision_svm <- svmConfMatrix$byClass['Pos Pred Value']
recall_svm <- svmConfMatrix$byClass['Sensitivity']
f1_svm <- 2 * (precision_svm * recall_svm) / (precision_svm + recall_svm)
precision_lda_Tweak1 <- ldaConfMatrixTweak1$byClass['Pos Pred Value']
recall_lda_Tweak1<- ldaConfMatrixTweak1$byClass['Sensitivity']
f1_lda_Tweak1 <- 2 * (precision_lda_Tweak1 * recall_lda_Tweak1) / (precision_lda_Tweak1 + recall_lda_Tweak1)
precision_svm_Tweak1 <- svmConfMatrixTweak1$byClass['Pos Pred Value']
recall_svm_Tweak1 <- svmConfMatrixTweak1$byClass['Sensitivity']
f1_svm_Tweak1 <- 2 * (precision_svm_Tweak1 * recall_svm_Tweak1) / (precision_svm_Tweak1 + recall_svm_Tweak1)
precision_svm_Tweak2 <- svmConfMatrixTweak2$byClass['Pos Pred Value']
recall_svm_Tweak2 <- svmConfMatrixTweak2$byClass['Sensitivity']
f1_svm_Tweak2 <- 2 * (precision_svm_Tweak2 * recall_svm_Tweak2) / (precision_svm_Tweak2 + recall_svm_Tweak2)
print(paste("GLM F1 Score:", f1_glm))
## [1] "GLM F1 Score: 0.607206142941524"
print(paste("LDA F1 Score:", f1_lda))
## [1] "LDA F1 Score: 0.591820368885325"
print(paste("SVM F1 Score:", f1_svm))
## [1] "SVM F1 Score: 0.779791323448655"
print(paste("LDA Tweak1 F1 Score:", f1_lda_Tweak1))
## [1] "LDA Tweak1 F1 Score: 0.000681663258350375"
print(paste("SVM Tweak1 F1 Score:", f1_svm_Tweak1))
## [1] "SVM Tweak1 F1 Score: 0.608448415922014"
print(paste("SVM Tweak2 F1 Score:", f1_svm_Tweak2))
## [1] "SVM Tweak2 F1 Score: 0.776057401812689"